#ifndef MATRIX_HPP_
#define MATRIX_HPP_

/*! \file matrix.hpp
 * \brief defines classes of matrices
 */ 

#include <stddef.h>
#include <stdexcept>
#include <vector>
#include <map>
#include <iterator>
#include "blinalg.hpp"
#include "usertypes.hpp"
using namespace std;

//*******************************

/*! \class Matrix
 * \author Guenther Of
 * \version 0.1
 * \date 16.04.2008
 * \brief Matrix is the base class of all kinds of matrices
*/
template<typename T> class Matrix
{
protected: 
  myindex nr_rows_; //!< number of rows
  myindex nr_cols_; //!< number of columns

public:
  virtual ~Matrix (){} //!< a
  myindex getNrRows () {return nr_rows_;} //!< \return number of rows
  myindex getNrCols () {return nr_cols_;} //!< \return number of columns
  virtual EVector<T> operator*(const EVector<T> &) const = 0; //!< matrix vector multiplcation is defined for deduced classes only 

};

//**********************************

/*! \class FullMatrix
 * \author Guenther Of
 * \version 0.1
 * \date 16.04.2008
 * \brief FullMatrix is a fully populated matrix
*/
template<typename T> class FullMatrix: public Matrix<T>
{
private:
  T** entry_; //!< memory for the matrix entries
public:
  FullMatrix (myindex n, myindex m, T t=0); //!< constructs n times m matrix, all entries set to t
  //FullMatrix (myindex n, myindex m, T** ptr); // n: #rows, m: #cols
  FullMatrix();
  FullMatrix (const FullMatrix &a); //!< copy constructor
  FullMatrix& operator=(const FullMatrix &a); //!< assignment operator 
  ~FullMatrix (); //!< destructor

  T* operator[] (myindex i) {return entry_[i];} //!< access to entry at row i and col j by [i][j]
  T& operator() (myindex i,myindex j) {return entry_[i][j];} //!< entry at row i and col j is (i,j)

  EVector<T> operator* (const EVector<T> &) const; //!< matrix vector multiplication 
  void setZero (); //!< all entries of the matrix are set to zero
  void print (); //!< matrix is printed

};

/*! \brief Constructor of FullMatrix, constructs a fully populated matrix of size n times m
 * \param[in] n number of rows
 * \param[in] m number of columns
 * \param[in] t intial value of all entries of the matrix (if dropped 0 is used)
 */
template<typename T> FullMatrix<T>::FullMatrix(myindex n, myindex m, T t) 
{
  this->nr_rows_ = n;
  this->nr_cols_ = m;

  entry_ = new T* [n];
  for (myindex i=0; i<n; i++)
  {
    entry_[i] = new T[m];

    for (myindex j=0; j<m; j++)
      entry_[i][j] = t;
  }
}

/*! \brief Copy constructor of FullMatrix, constructs a fully populated matrix as copy of a
 * \param[in] a matrix to be copied
 */
template<typename T> FullMatrix<T>::FullMatrix(const FullMatrix &a)
{
  this->nr_rows_ = a.nr_rows_;
  this->nr_cols_ = a.nr_cols_;

  entry_ = new T* [this->nr_rows_];
  for (myindex i=0; i<this->nr_rows_; i++)
  {
    entry_[i] = new T[this->nr_cols_];

    for (myindex j=0; j<this->nr_cols_; j++)
      entry_[i][j] = a.entry_[i][j];
  }
}
template<typename T> FullMatrix<T>::FullMatrix()
{
	FullMatrix(0, 0, 0);
}
/*! \brief Assingment operator of FullMatrix, 
 * \param[in] a matrix to be assigned
 * \return copy of matrix a
 * \exception std::out_of_range matrix dimensions do not match
 */
template<typename T> 
FullMatrix<T> & FullMatrix<T>::operator=(const FullMatrix &a)
{
  if (this != &a)
  {
    if ((this->nr_rows_ != a.nr_rows_) || (this->nr_cols_ != a.nr_cols_))
      throw std::out_of_range ("matrix sizes do not match in FullMatrix<T>::operator=()");

    for (myindex i=0; i<this->nr_rows_; i++)
      for (myindex j=0; j<this->nr_cols_; j++)
	entry_[i][j] = a.entry_[i][j];
  }
  return *this;
}

/*! \brief Destructor of FullMatrix
 */
template<typename T> FullMatrix<T>::~FullMatrix()
{
  for (myindex i=0; i<this->nr_rows_; i++)
    delete [] entry_[i];

  delete[] entry_;
}

/*! \brief matrix times vector multiplication z = A v
 * \param[in] v vector
 * \return result of matrix times vector product A v
 * \exception std::out_of_range dimensions of the matrix and of the vectors do not match
 */
template<typename T>
EVector<T> FullMatrix<T>::operator* (const EVector<T> &v) const
{
  if (this->nr_cols_ != v.size())
    throw std::out_of_range ("matrix and vector sizes do not match in FullMatrix<T>::operator*");

  EVector<T> res(this->nr_rows_);

  for (myindex i=0; i<this->nr_rows_; i++)
  {
    res[i] = 0;
    for (myindex j=0; j<this->nr_cols_; j++)
      res[i] += entry_[i][j]*v[j];
  }

  return res;
}

/*! \brief sets all entires of the matrix to zero.
 */
template<typename T>
void FullMatrix<T>::setZero ()
{
  for (myindex i=0; i<this->nr_rows_; i++)
  {
    for (myindex j=0; j<this->nr_cols_; j++)
      entry_[i][j] = 0;
  }
}

/*! \brief prints the matrix to cout
 */
template<typename T>
void FullMatrix<T>::print ()
{
  for (myindex i=0; i<this->nr_rows_; i++)
  {
    for (myindex j=0; j<this->nr_cols_; j++)
//      printf ("%1.17e ", entry_[i][j]);
      std::cout << entry_[i][j] << " ";

    printf ("\n");
      //std::cout << ";" << std::endl;
  }
}

//*************************

/*! \class SymmetricFullMatrix
 * \author Guenther Of
 * \version 0.1
 * \date 16.04.2008
 * \brief SymmetricFullMatrix is a symmetric fully populated matrix, stored as lower triangular matrix
*/
template<typename T> class SymmetricFullMatrix: public Matrix<T>
{
private:
  T** entry_; //!< memory for the matrix entries
public:
  SymmetricFullMatrix (myindex n, T t=0);  //!< constructs symmetric n times m matrix, all entries set to t
  SymmetricFullMatrix (const SymmetricFullMatrix &a); //!< copy constructor
  SymmetricFullMatrix& operator=(const SymmetricFullMatrix &a); //!< assignment operator
  ~SymmetricFullMatrix (); //!< destructor
 
  T* operator[] (myindex i); //!< access to entry at row i and col j by [i][j]
  T& operator() (myindex i,myindex j); //!< entry at row i and col j is (i,j)

  EVector<T> operator* (const EVector<T> &) const; //!< matrix vector multiplication
  void setZero (); //!< all entries of the matrix are set to zero
  void print (); //!< matrix is printed
};

/*! \brief Constructor of SymmetricFullMatrix, constructs a symmetric fully populated matrix 
 * of size n times m
 * \param[in] n number of rows and columns
 * \param[in] t intial value of all entries of the matrix (if dropped 0 is used)
 */
template<typename T> SymmetricFullMatrix<T>::SymmetricFullMatrix(myindex n, T t) 
{
  this->nr_rows_ = n;
  this->nr_cols_ = n;

  entry_ = new T* [n];
  for (myindex i=0; i<n; i++)
  {
    entry_[i] = new T[i+1];

    for (myindex j=0; j<i+1; j++)
      entry_[i][j] = t;
  }
}

/*! \brief Copy constructor of SymmetricFullMatrix, constructs a symmetric fully populated matrix 
 * as copy of the matrix a
 * \param[in] a matrix to be copied
 */
template<typename T> SymmetricFullMatrix<T>::SymmetricFullMatrix(const SymmetricFullMatrix &a)
{
  this->nr_rows_ = a.nr_rows_;
  this->nr_cols_ = a.nr_cols_;

  entry_ = new T* [this->nr_rows_];
  for (myindex i=0; i<this->nr_rows_; i++)
  {
    entry_[i] = new T[i+1];

    for (myindex j=0; j<i+1; j++)
      entry_[i][j] = a.entry_[i][j];
  }
}

/*! \brief Assingment operator of SymmetricFullMatrix, 
 * \param[in] a matrix to be assigned
 * \return copy of matrix a
 * \exception std::out_of_range matrix dimensions do not match
 */
template<typename T> 
SymmetricFullMatrix<T> & SymmetricFullMatrix<T>::operator=(const SymmetricFullMatrix &a)
{
  if (this != &a)
  {
    if ((this->nr_rows_ != a.nr_rows_) || (this->nr_cols_ != a.nr_cols_))
      throw std::out_of_range ("matrix sizes do not match in FullMatrix<T>::operator=()");

    for (myindex i=0; i<this->nr_rows_; i++)
      for (myindex j=0; j<i+1; j++)
	entry_[i][j] = a.entry_[i][j];
  }
  return *this;
}

/*! \brief Destructor of SymmetricFullMatrix
 */
template<typename T> SymmetricFullMatrix<T>::~SymmetricFullMatrix()
{
  for (myindex i=0; i<this->nr_rows_; i++)
    delete [] entry_[i];

  delete[] entry_;
}

/*! \brief provides access to entry at row i and col j by [i][j]
 * \warning Access the lower traingular matrix only. 
 * If you access to th upper part, anything (crash) can happen
 */
template<typename T> 
T* SymmetricFullMatrix<T>::operator[] (myindex i) {return entry_[i];}

/*! \brief provides access to entry  entry at row i and col j is (i,j).
 * You can access any entry of the full matrix
 */
template<typename T> 
T& SymmetricFullMatrix<T>::operator() (myindex i,myindex j) 
{
  if (j>i)
    return entry_[j][i];
  else
    return entry_[i][j];
}

/*! \brief matrix times vector multiplication z = A v of a symmetric matrix 
 * \param[in] v vector
 * \return result of matrix times vector product A v
 * \exception std::out_of_range dimensions of the matrix and of the vectors do not match
 */
template<typename T>
EVector<T> SymmetricFullMatrix<T>::operator* (const EVector<T> &v) const
{
  if (this->nr_cols_ != v.size())
    throw std::out_of_range ("matrix and vector sizes do not match in FullMatrix<T>::operator*");

  EVector<T> res(this->nr_rows_);

  for (myindex i=0; i<this->nr_rows_; i++)
  {
    res[i] = 0;
    for (myindex j=0; j<i; j++)
    {
      res[i] += entry_[i][j]*v[j];
      res[j] += entry_[i][j]*v[i];
    }
    res[i] += entry_[i][i]*v[i];
  }

  return res;
}

/*! \brief sets all entires of the matrix to zero.
 */
template<typename T> 
void SymmetricFullMatrix<T>::setZero()
{
  
  for (myindex i=0; i<this->nr_rows_; i++)
    for (myindex j=0; j<i+1; j++)
      entry_[i][j] = 0;
}

/*! \brief prints the matrix to cout
 */
template<typename T>
void SymmetricFullMatrix<T>::print ()
{
  for (myindex i=0; i<this->nr_rows_; i++)
  {
    for (myindex j=0; j<i; j++)
      std::cout << entry_[i][j] << " ";
    for (myindex j=i; j<this->nr_cols_; j++)
      std::cout << entry_[j][i] << " ";

    std::cout << std::endl;
  }

}

//******************************
/*! \class SparseMatrix
 * \author Guenther Of
 * \version 0.1
 * \date 16.04.2008
 * \brief SparseMatrix is a sparse matrix
*/
template<typename T> class SparseMatrix: public Matrix<T>
{
private:
  std::map<myindex,T>* entry_; //!< memory for the matrix entries
 
public:
  SparseMatrix (myindex n, myindex m); ///!< constructs n times m matrix, all entries set to t
  //SparseMatrix (myindex n, T** ptr); // n: #rows, m: #cols
  SparseMatrix (const SparseMatrix &a); //!< copy constructor
  SparseMatrix& operator=(const SparseMatrix &a); //!< assignment operator
  ~SparseMatrix (); //!< destructor
 
  std::map<myindex,T>& operator[] (myindex i); //!< access to entry at row i and col j by [i][j]
  T& operator() (myindex i,myindex j); //!< entry at row i and col j is (i,j)

  EVector<T> operator* (const EVector<T> &) const;  //!< matrix vector multiplication
  void setZero (); //!< all entries of the matrix are set to zero
  void print (); //!< print matrix to cout
};


/*! \brief Constructor of SparseMatrix, constructs a sparse matrix of size n times m
 * \param[in] n number of rows
 * \param[in] m number of columns
 */
template<typename T> SparseMatrix<T>::SparseMatrix(myindex n, myindex m) 
{
  this->nr_rows_ = n;
  this->nr_cols_ = m;

  entry_ = new std::map<myindex,T> [n];
}

/*! \brief Destructor of SparseMatrix
 */
template<typename T> SparseMatrix<T>::~SparseMatrix()
{
  delete[] entry_;
}

/*! \brief Copy constructor of SparseMatrix, constructs a sparse matrix as copy of matrix a
 * \param[in] a matrix to be copied
 */
template<typename T> SparseMatrix<T>::SparseMatrix(const SparseMatrix &a)
{
  this->nr_rows_ = a.nr_rows_;
  this->nr_cols_ = a.nr_cols_;

  entry_ = new std::map<myindex,T> [this->nr_rows_];
  for (myindex i=0; i<this->nr_rows_; i++)
    entry_[i] = a.entry_[i];
}

/*! \brief Assingment operator of SparseMatrix, 
 * \param[in] a matrix to be assigned
 * \return copy of matrix a
 * \exception std::out_of_range matrix dimensions do not match
 */
template<typename T> 
SparseMatrix<T> & SparseMatrix<T>::operator=(const SparseMatrix &a)
{
  if (this != &a)
  {
    if ((this->nr_rows_ != a.nr_rows_) || (this->nr_cols_ != a.nr_cols_))
      throw std::out_of_range ("matrix sizes do not match in FullMatrix<T>::operator=()");

    for (myindex i=0; i<this->nr_rows_; i++)
      entry_[i] = a.entry_[i];
  }
  return *this;
}

/*! \brief provides access to entry at row i and col j by [i][j]
 * If the entries doesnot exist, it will be created.
 */
template<typename T> 
std::map<myindex,T>& SparseMatrix<T>::operator[] (myindex i) {return entry_[i];}

/*! \brief provides access to entry  entry at row i and col j is (i,j).
 * If the entries doesnot exist, it will be created.
 */
template<typename T> 
T& SparseMatrix<T>::operator() (myindex i,myindex j) 
{return entry_[i][j];}

/*! \brief matrix times vector multiplication z = A v for a sparse matrix A
 * \param[in] v vector
 * \return result of matrix times vector product A v
 * \exception std::out_of_range dimensions of the matrix and of the vectors do not match
 */
template<typename T>
EVector<T> SparseMatrix<T>::operator* (const EVector<T> &v) const
{
  if (this->nr_cols_ != v.size())
    throw std::out_of_range ("matrix and vector sizes do not match in FullMatrix<T>::operator*");

  std::map<myindex,double>::iterator p; //todo: does not work with T
  EVector<T> res(this->nr_rows_);

  for (myindex i=0; i<this->nr_rows_; i++)
  {
    res[i] = 0;
    for (p = entry_[i].begin(); p != entry_[i].end(); p++)
      res[i] += p->second*v[p->first];
  }

  return res;
}

/*! \brief sets all entires of the matrix to zero.
 */
template<typename T>
void SparseMatrix<T>::setZero ()
{
  for (myindex i=0; i<this->nr_rows_; i++)
    entry_[i].clear();
}

/*! \brief prints the matrix to cout
 */
template<typename T>
void SparseMatrix<T>::print ()
{
  std::map<myindex,double>::iterator p; //todo: does not work with T

  for (myindex i=0; i<this->nr_rows_; i++)
  {
    std::cout << "row " << i <<": ";
    for (p = entry_[i].begin(); p != entry_[i].end(); p++)
      std::cout << p->first << ": " << p->second << ", ";

    std::cout << std::endl;
  }
}

#endif // MATRIX_HPP_
